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ABSTRACT 


A  laminar  wake  model  is  studied  which  assumes  that  the  boundary 
layer  equations  are  valid  in  a  shear  layer  on  the  axis  of  symmetry,  and 
that  the  local  pressure  is  determined  by  a  displacement  interaction  with 
the  inviscid  outer  flow.  The  finite  difference  method  used  places  no 
restrictions  on  the  shape  of  the  velocity  profile.  The  qualitative 
features  of  the  solution  are  the  same  as  those  found  by  previous  investi¬ 
gators  using  integral  methods,  in  that  an  otherwise  arbitrary  initial 
conation  (the  local  free  stream  Mach  number)  is  determined  by  the  fact 
that  only  a  particular  value  results  in  a  physically  meaningful  solu¬ 
tion.  Two  methods  of  finding  the  particular  solution  are  investigated. 


This  research  is  a  part  of  PROJECT  DEFENDER,  sponsored  by  the  Advanced 
Research  Projects  Agency,  Department  of  Defense,  under  ARPA  Order  No. 

AF04(694)1-570Ored  ^  ^  BalUstic  Syscems  Division  under  Contract 


i 

EOS  RN-28,  10-65 


i 


TABLE  OF  CONTENTS 


INTRODUCTION 

EQUATIONS 

METHOD  OF  SOLUTION 
RESULTS 

APPENDIX  I  -  INVERSE  METHOD 
FIGURES 


1 

5 

10 

13 

16 

19 


EOS  RN-28,  10-65 


TirOS  PACE  WAS  BIAHK,  TH3HSFCS3  NOT  FI  IKED, 


ILLUSTRATIONS 


1.  Base  flow  region  19 

2«  du*/dS*  for  M  close  to  the  eigenvalue  20 

c  eo 

3*  dlnM^/dSA  for  M  close  to  the  eigenvalue  20 

4.  Comparison  between  profile  shapes  from  the  nor-similar  21 

calculation  and  those  from  the  Stewar tson-Kennedy 
similar  solution 


5. 

Rate  of  convergence 

to  the  eigenvalue 

21 

6. 

Convergence  of  Inverse  method,  20  iterations 

21 

7. 

Comparison  of  direct 

and  inverse  methods,  du*/dS 

22 

8. 

Comparison  of  direct 

and  inverse  methods,  dlnMe/dS* 

22 

EOS  RN-28,  10-65  v 

ERSVI0D3  PAGE  WAS  HLAJ2C,  THSSSTORK  NOT  FUMED* 


* 

m 


ACKNOWLEDGMENT 


The  author  wishes  to  express  his  appreciation  to  M.R .  Denison 
for  his  encouragement  and  for  many  useful  discussions  during  the  c 
of  this  work. 


se 


EOS  RN-28,  10-65 


vii 


A  Finite  Difference  Method  for  Solving  the  Laminar 
Boundary  Layer  Equations  Near  a  Wake  Neck 


INTRODUCTION 

The  qualitative  features  of  the  base  flow  region  of  supersonic  blunt- 
based  bodies,  shown  schematically  in  Fig.  L  have  been  recognized  for 
many  years.  Much  of  the  early  work  was  aimed  at  determining  the  base 
pressure,  and  it  was  generally  agreed,  long  before  any  quantitative  theory 
was  available,  that  an  interaction  between  the  shear  layer  and  the  outer 
inviscid  flow  was  important  in  determining  its  value.  The  Crocco-Lees^ 
mixing  theory  put  this  on  a  somewhat  more  quantitative  basis.  The  model 
on  which  this  theory  is  based  assumes  that  the  boundary  layer  equations 
are  valid  in  the  shear  layer  and  that  the  local  pressure  Is  determined  by 
a  displacement  interaction  with  the  inviscid  outer  flow.  The  effects  of 
the  shear  layer  profiles  are  lumped  into  variables  which  must  be  deter¬ 
mined  independently,  so  that  the  theory  is  incomplete.  However,  the 
qualitative  results  show  what  behavior  to  expect  from  a  more  rigorous 
solution  of  the  same  model.  They  found  that  there  were  an  infinite 
number  of  solutions  to  the  equations,  corresponding  to  all  possible 
choices  of  the  base  pressure.  However,  the  only  one  with  physical  mean¬ 
ing  was  the  singular  solution  passing  through  a  saddle-point  type  of 
singularity.  The  most  familiar  example  of  this  type  of  behavior  is  in 
the  calculation  of  flow  through  a  one  dimensional  nozzle.  In  that  case, 
the  governing  equation  is 

dM  i  dA  MQ  +  ^ 

dX  '  '  A  d*  (1  -  M2) 


which  can  of  course  be  integrated  in  closed  form.  However,  if  one 

1  dA 

chooses  to  consider  ~  ^  as  a  specified  function  of  x,  one  can  solve 
the  first  ovder  differential  equation  for  any  initial  value  of  M. 


1.  L.  Crocco  and  L.  Lees,  MA  Mixing  Theory  for  the  Interaction  Between 
Dissipative  Flows  and  Nearly  Isentropic  Streams, "  J.  Aer.  Sci.  19, 
649-676  (1952).  “ 
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If  the  required  solution  is  to  have  M  <  1  in  the  converging  section  and 

dA 

M  y  1  in  the  diverging  section,  then  M  must  pass  through  1  and  — 

through  zero  during  the  course  of  the  calculation-  If  the  initial  choice 

,  dA  dM 

oi  M  is  too  small,  —  will  pass  through  zero  first,  after  which 

becomes  negative  and  the  flow  is  subsonic  everywhere  in  the  nozzle. 

If  the  initial  choice  is  too  large  M  -  1  first  and  the  calculation  ter¬ 
minates  at  M  =  1,  a  solution  with  no  physical  meaning.  The  only  solu¬ 
tion  which  becomes  supersonic  in  the  diverging  section  corresponds  to 
“  *0  and  M  -  1  simultaneously. 


This  analogy  is  useful  in  obtaining  an  understanding  of  the  nature 
of  the  wake  critical  point,  but  it  should  not  be  carried  too  far.  For 
example,  the*  pressure  increases  as  the  flow  approaches  the  wake  critical 
point,  indicating  that  the  acceleration  of  the  low  speed  portion  of  the 
flow  is  produced  by  viscous  f orxes ,  which  are  not  present  in  the  one 
dimensional  nozzle  model.  Relating  the  wake  critical  point  to  an  aver¬ 
age  Mach  number  of  unity  is  therefore  not  justified. 


Recently,  Reeves  and  Lees  investigated  this  wake  model  in  more  detail, 

using  a  two  moment  integral  method  to  determine  profile  shape  and  wake 

3 

thickness.  To  insure  the  proper  qualitative  features,  the  Stewartson 
famiLy  of  solutions  to  the  Falkner-Skan  eq  'Lion,  f1  (^-7^),  was  used  to 
describe  the  profile  shape,  ~  (Jj— ,y)  with  ij^  and  y  determined  by  the 
zero  and  firsc  moments  of  the  momentum  equation.  Since  a  stagnation 
point  on  the  centerline  occurs  for  only  a  single  value  of  y,  the  shape 
of  the  profile  at  the  wake  stagnation  point  is  arbitrarily  predeter¬ 
mined  by  the  choice  of  the  Stewartson  family  and  only  the  wake  thickness 
at  this  point  is  calculated. 


2.  B,  L.  Reeves  and  L.  Lees,  "Theory  of  the  Laminar  Near  Wake  of 
Blunt  Bodies  in  Hypersonic  Flow”  presented  at  AIAA  2nd  Aerospace 
Sciences  Meeting,  New  York,  January  1965,  AIAA  Paper  No.  6 5 •  •  5 2 . 

3.  K,  Stewartson,  further  Solutions  of  the  Falkner  Skan  Equation," 
Proc .  Comb „  Phil.  Soc.  50,  454-465  (1954). 


EOS  RN-28 ,  10-65 


1 


J 


1 


L 


2 


Webb,  Golik,  Vogenitz  and  Lees  investigated  the  model  using  the 
polynomial  profiles 


u_ 

u. 


and 


u 

uc 


+  1  " 


in  which  a  two  moment  method  was  used  to  find 


and 


He  e,  again 


the  profile  shape  at  the  wake  stagnation  point  is  predetermined  by  the 

choice  of  profile  families.  They  then  went  one  step  further  and  used 

T 

the  Stewartson  family  f 1  (*yp,  p)  to  describe  profile  shape,  using  the 
following  substitutions: 

f  -  (u-uc)/(u&-uc) 


P  -  Y 

u 

c 

where  —  is  now  no  longer  related  to  y.  The  third  parameter  was  found 
using  an  additional  moment  of  the  momentum  equ  tion. 

The  biggest  drawback  of  the  integral  methods  is  that  in  order  to 
get  accurate  results,  the  family  of  profiles  chosen  mus  be  capable  of 
a  reasonably  good  representation  of  the  actual  profiles  over  the  entire 
range  of  calculation.  Since  the  actual  profiles  are  not  generally  known 
in  advance,  this  requires  the  exercise  of  a  considerable  degree  of  judg¬ 
ment  and  intuition  on  the  part  of  the  investigator.  In  this  particular 
problem,  there  is  the  additional  obstacle  of  a  gap  between  the  body 
boundary  layer  and  the  neighborhood  of  the  rear  stagnation  point,  in 
which  the  boundary  layer  equations  alone  are  clearly  inadequate  in 
describing  the  flow  between  the  centerline  and  the  outer  inviscid  flow. 


4.  W.H.  Webb,  R.J.  Golik,  F.W.  Vogenitz  and  Lester  Lees,  ,fA  Multimoment 
Integral  Theory  of  the  Laminar  Supersonic  Near  Wake,”  Proc . ,  1965 
Heat  Transfer  and  Fluid  Mechanics  Institute,  p.  168,  Stanford  U.  Pr~ss, 
Stanford,  Calif-,  1965. 
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Since  a  different  method  of  calculation  is  required  iu  this  region,  it 
is  not  likely  that  the  profiles  approaching  the  rear  stagnation  point 
will  fit  smoothly  intc  a  precnosen  family,  and  it  is  not  possible  to 
an  integral  method  calculation  on  an  arbitrary  initial  profile. 

Hits  paper  describes  a  finite  difference  method  for  solving  the 
boundary  layer  equations  with  the  local  pressure  determined  by  a  displace¬ 
ment  interactiou  with  the  outside  inviscid  flow.  There  are  no  constraints 
on  profile  shape,  and  the  calculation  begins  with  an  arbitrary  profile 
at  the  wake  stagnation  point  and  proceeds  downstream.  The  solutions  r.re 
found  to  fall  into  two  families,  depending  on  whether  the  initial  loer 1 
inviscid  Mach  number  is  smaller  or  larger  than  the  value  corresponding 
to  a  singular  solution.  If  the  Mach  number  is  smaller  than  this  eigen¬ 
value,  the  centerline  velocity  increases  ror  only  a  short  distance 

downstream,  then  decreases  again  until  another  stagnation  point  is 
reached,  the  flow  past  this  point  being  in  the  direction  of  the  body. 

If  the  initial  local  inviscid  Mach  number  is  larger  than  the  eigenvalue, 
the  shear  layer  thickness  decreases  at  an  accelerating  rate,  until  the 
inviscid  flow  intersects  the  axis  as  the  shear  layer  disappears.  Both 
of  these  types  of  behavior  are  valid  solutions  to  the  equations,  but 
only  the  singular  solution  has  the  appropriate  downstream  behavior, 
that  is,  the  wake  velocity  approaches  the  local  freestream  value. 
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EQUATIONS 

The  coordinate  system  as  shown  in  Fig.  1  employs  distances  measured 
along  the  centerline  and  normal  to  it.  The  origin  is  at  the  base  flow 
rear  stagnation  point.  The  calculations  will  be  limited  here  to  two 
dimensional,  adiabatic,  ideal  gas,  flow.  The  continuity  and  momentum 
equations 

5pu  Spv  A  >  k 

sr +  §r  ■ 0  CD 


PU 


bu 

bx 


+  pv 


du 


b  bu 
by  by 


dP 

dx 


(2) 


can  be  transformed 


u* 


bu* 

bs* 


4-  v* 


into  the  form 

bxx *  52u* 

dY*  ~  Sy*?  * 


dP 

dS*  * 


1 


p  u  * 
e  e 


0) 


where  u*  =  u/u^,  e  denoting  conditions  at  shear  layer  edge. 


o 


S* 


1 


Cp  u  \x  dx  where  C  -  — - 
e  e  e  p  }i 

He  e 


(Chapman -Rube sin  parameter, 
assumed  to  be  a  function  of 
x  only) 
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=  scale  factor  making  equation  dimensionless 


For  ideal  gas  conditions,  the  right  side  of  Eq.  (3)  can  be  written 

dM 


dP 

dS*  p 


P 

i  e  __ 

IT2  p 

i  e  L 


H  dS* 
e 


e  f  i  2\ 

(!  -  u*  ) 


(4) 


It  is  convenient,  now,  to  further  transform  Eqp  (3)  into  a  less  familiar 
form,  with  dependent  variable 

^2 

(5) 


and  independent  variables  S*and  u,  where 


u*  -  U* 

—  Q 

u  =  - - — c  denoting  conditions  at  the  centerline 

i  - 

c 


(6) 


M 


dG 

ds* 


_(1  -  u*)2}  dG2 

c 


c^G  ,  |  ~u*(l  -  i 

2  a,:2  La  -  «*) 


„  du* 

-  11*)  c 


gg .  i-  S§  (i-u*2) 


*)  dS*  2(l-u*)  Su  Me  dS* 


1  1 

,  dG 

Q 

dM  H 

_ e_  „* 

-  u*) 

du 

M 

dS* 

C  ^ 

L  e 

(7) 


The  boundary  conditions  appropriate  to  a  wake  are  G(o)  *  G(l)  «  0# 

The  fixed  location  of  the  boundaries  in  the  u  coordinate  greatly 

simplify  the  application  of  the  finite  difference  method,  and  the  use 
£u* 

of  G  instead  of  as  dependent  variable  makes  it  more  convenient  to 

satisfy  the  differential  equation  at  the  centerline,  where  Eq*  (7) 
reduces  to 

*  Hi  1  I  Sg\  ,  _1  Hi  n  *2^ 

Uc  dS*  "  2(1  -  u*)  1  du  !  +  M  dS*  "  uc  * 

'  c '  \  /c  e 
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Here  since  i  *  f  j  y  ,  the  advantage  in  using  G  as  the 

\  UA  v  /c  V  u  4  /\*\ 

dependent  variable  becomes  clear  when  it  is  realized  that,  while  |  r-r  I 

f  \  (  \  vY  M 

is  zero,  (  *“S  ^n^^n^-te  so  that  f-4jj  generally  has  a  finite 

value , 

Equation  (7)  is  in  a  modified  form  of  the  Crocco  coordinate  system 
and,  as  is  well  known,  solutions  to  the  boundary  layer  equations  in  these 
coordinates  are  not  unique,  in  that  a  single  solution  in  the  Crocco 
coordinates  contains  an  infinite  number  of  distinct  solutions  in  the 
physical  coordinates.  This  occurs  because  the  usual  third  boundary 
condition  (the  location  of  a  reference  streamline)  of  the  momentum 
equation  in  physical  coordinates  cannot  be  imposed  in  the  Crocco  coordin¬ 
ates,  and  the  solutions  corresponding  to  all  possible  values  of  this 
boundary  condition  map  into  a  single  Crocco  solution.  We  wish,  here,  to 
remove  the  ambiguity  of  the  solution  of  Eq.  7  by  imposing  the  condition 
thaw  the  centerline  (u  =  0)  is  a  streamline.  Of  the  many  possible  ways 
of  accomplishing  this,  we  choose  to  use  an  integral  of  the  momentum  equa¬ 
tion,  with  (^c\  set  equal  to  zero  at  Y  =  0: 


a  -  u*>  \  u*  ci .  u*)  ^=1  =  - 

L  \  i  g  J 


-  u*)\  (1 


U*2)^ 


An  additional  relation  between 


is  required  i“ 


order  to  complete  the  specification  of  the  problem.  This  is  obtained 
from  an  integrated  form  of  the  continuity  equation: 


d  in  p^u^ 


tan  9 


where 


EOS  RN-28 ,  10-65 


The  logarithmic  term  in  Eq.  (9)  is  not  bounded  with  increasing 
5,  and  since  5  is  arbitrary,  the  inclusion  of  this  term  is  difficult  to 
justify  within  the  framework  of  the  boundary  layer  theory,  and  it  will 
be  neglected  in  these  calculations. 

Equation  (9),  when  transformed  into  the  coordinate  system  of 


Eq.  (3),  becomes 


p  u  6* 

*00  CD 


M  /  1  +  ¥ 

Me  \  1  +  ^  M, 


tan  0 


where  it  has  been  assumed  for  simplicity  that  |i  ~  T  '  .  Subscript 
refers  to  a  reference  point  in  the  flow,  tiken  here  to  be  at  0  *  0* 

'FT  r— 

The  term  —  is  like  \ Re „  If  the  calculation  were  linked  to  a  given 

f^co 

body,  it  would  be  convenient  to  divide  this  into  a  conventional 
Reynolds  number  and  a  factor  dependent  primarily  on  body  shape,  but 
for  the  present  we  will  define 

sb 

Re  =  — >j —  fl 


A  transformation  of  the  displacement  thickness  definition,  Eq.  (10) 
into  the  coordinate  system  of  Eq.  (3)  gives: 


(13) 

+ 


where,  it  should  be  noted,  both  of  the  shear  layer  integral  functions 
are  bounded . 
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A  final  relation  is  required  between  the  local  flow  angle  9  and 

local  inviscid  Mach  number  M  .  This  could  be  obtained  by  means  of  a 

e 

detailed  calculation  of  the  inviscid  flow  field,  but  for  the  purposes 
of  these  calculations,  we  will  assume  that  this  relation  is  given  by 
the  Prandtl-Meyer  function: 

Q  -  -  js  (14) 

<x 


The  equations  presented  in  this  section,  along  with  the  appropriate 
initial  conditions,  completely  specify  the  problem. 

The  momentum  boundary  layer  equation  (7)  requires  an  initial  G  vs  u 


profile,  and  auxilliary  equations  to  evaluate 

Equations  (8),  (11),  (13)  and  (14)  relate 

d/p_um&*^ 


-  Me 


du* 
_ c 

dS*’ 


.  dM 
l  a 

M  dS*’ 


dS* 


u*,  M  , 
c>  e> 


p  u  5* 

Kco  CO 


fs: 


dMe  duQ 

dS*  and  dS* 


and  0 


If  0  is  eliminated  between  Eqs.  (14)  and  (11)  and  equation  (13)  is 
differentiated  the  three  remaining  equations  could  be  thought  of  as 

P<X)^CO^^ 

simultaneous  ordinary  differential  equations  for  u*.  M  and  — - f 

c  e  rh 

where  only  two  initial  conditions  are  arbitrary  because  Eq#  (13)  pro¬ 
vides  a  relation  between  these  and  the  third.  Since  one  would  normally 


know  u*  if  he  knows  the  initial  G  vs.  u  profile,  we.will  consider 
c  f 

(u*) initial  to  be  and  then  arbitrarily  1st  ( - 


be 


'initial 

specified  by  Eq.  (13),  leaving  (M^)  ^n^t£a]  as  Yet  undetermined.  MeQ 

will  then  be  the  eigenvalue  and  will  be  chosen  so  that  the  solution  far 
downstream  has  the  appropriate  behavior. 
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METHOD  OF  SOLUTION 


The  momentum  equation  (7)  can  be  solved  readily  by  the  use  of  an 
implicit  finite  difference  method,  with  the  accuracy  limited  only  by 
the  available  computing  time  and  computer  memory  size.  Initial  condi¬ 
tions  are  unrestricted  except  for  the  requirement  that  they  satisfy  the 
boundary  conditions  and  not  contain  discontinuities.  Exceptions  to 
these  requirements  can  be  successfully  treated  by  using  local  analytical 
solutions,  but  we  do  not  expect  any  such  cases  to  arise  in  the  current 
problem.  The  method  uses  the  following  finite  difference  approximations 
for  the  linear  (not  bracketed)  terms  in  Eq.  (7)  : 


9G 

c)S 


(Gi,  j+i’ 


Gij> /AS 


(15) 


CG,L+l,j+l+  Gi+l,j' 


<Gi,j+l+Gi,j>  +Gi-i,j+l+Gi-l,j>/2AG 


Sg 

du 


^Gi+lJ+l+  Gi+l,j"  Gi-1, j+1 


C1-1J>  '  4AG 


(17) 


G 


(18) 


where  i  refers  to  the  u  coordinate 
j  refers  to  the  S  coordinate 

The  bracketed  terms  are  evaluated  at  i,j  for  an  initial  iteration,  the.i 
averaged  between  i,j  and  i,j*fl  for  succeeding  iterations  until  they  con¬ 
verge  to  a  limiting  value.  If  N  interior  meshpoints  in  u  are  used,  the 
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finite  difference  approximation  to  Equation  (7)  results  in  N  simultan¬ 
eous  linear  algebraic  equations  for  the  unknown  G.  ....  The  equation 

i ,  j+1 

at  each  point  refers  to  Gi+1  ^  and  G^  +1  and  therefore  the  equa¬ 
tions  are  coupled  in  such  a  way  that  the  matrix  of  coefficients  contains 
3  bands  consisting  of  the  diagonal,  and  the  bands  immediately  above  and 
below  the  diagonal.  Such  a  set  of  equations  is  very  quickly  and  easily 
solved  by  Gaussian  elimination. 

The  auxiliary  equations  are  solved  simultaneously  with  the  momentum 
equation  in  the  following  sequence. 


1.  Guess 


ss*h+2 


2.  Calculate  u* 
c 


3.  Solve  the  momentum  equation ,dgvaluating  bracketed  terras  at  j 
and  using  the  previous  g-  “•  for  the  first  iteration. 

fl  dMe\  6 

4.  Calculate  /—  ~J  ,  from  Eq.  (8). 


5.  Repeat  three  and  four  using  averaged  values  until  converged. 


dM 

6.  Calculate  Me  by  finite  difference  integration  of  — $■ 

/ ...  MQ  ds 


7 .  Calculate 


o  u_6* 

CO  00 


from  Eq.  (13). 


Jb  /j+1 


A  f  p  u  6*\ 

8.  Evaluate  — ..[  — — ” —  ! 

s  \  is-/  , 

'  ’  b  'j+1/ 2 


from  Eq.  (11) 


9.  Compare  the  finite  difference  derivative 

7p»u»6*\  (pm*x 

Irav  lirir 

with  that  evaluated  in  Step  8. 
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'du*\ 


10.  Use  a  systematic  iteration  procedure  to  guess  a  new  i 


\  °  /i+1/2 

which  will  give  better  agreement  than  in  step  9  and  repeat' 
from  Step  1  to  convergence. 


(dut\ 

The  modified  Newtons  method  used  for  iterating  on  l  dS*  ] 
converges  to  5  figures  in  about  4  iterations. 


/J+l/2 


generally 


The  principal  drawback  of  the  method  is  due  to  Step  9,  where  a  finite 
difference  derivative  is  required.  In  order  for  this  difference  to  main¬ 
tain  significance,  the  profiles  which  determine  5*  must  be  considerably 
more  accurate  then  one  would  normally  require  and  the  calculation  there¬ 
fore  more  time-consuming. 
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RESULTS 


Typica1  results  of  this  method  of  calculation  are  shown  in  Figs, 

(2)  and  (3),  The  initial  profile  used  was,  for  convenience,  one  of 

the  Stewarcson  family  tabulated  by  Kennedy,5  where  f"  ^ »P)~  \G~(u, 1988) 

with  Ro-l(f  and  Mg=6 .  Figure  (2)  shows  the  behavior  of  -r~~  corresponding  to 
different  assumed  eigenvalues.  When  Meo  is  too  large,  —begins  to 
increase  rapidly  atsogg  point  in  the  calculation.  This  is  accompanied 
by  an  increase  in  j-j-  toward  infinity  from  its  initially  negative 
value,  as  shown  in  Fig.  (3).  The  centerline  velocity  is  therefore 
increasing  rapidly  at  the  same  time  that  the  inviscid  flow  just  outside 
the  shear  layer  is  turning  rapidly  in  toward  the  axis  of  symmetry.  In 
the  limit,  the  shear  layer  disappears  entirely  as  the  inviscid  flow 
intersects  the  axis.  This  limit  is  reached  at  a  finite  distance  down¬ 
stream  and  constitutes  a  barrier  past  which  the  calculgtion  cannot  pro¬ 
ceed.  When  the  assumed  eigenvalue  is  too  small,  becomes  nega¬ 

tive  at  some  point  in  the  calculation,  after  which  the  centerline  velocity 
ratio  decreases  until  a  second  stagnation  point  is  reached,  the  flow  on 

the  centerline  past  this  point  having  its  origin  downstream  rather  than 
upstream. 

Only  the  eigen-  solution  belongs  to  neither  of  these  two  families 
of  non-physical  solutions.  The  source  of  this  singular  behavior  can  be 
more  easily  seen  by  differentiating  equatior  (13),  and  combining  with 
Eq.  (11)  to  get  an  equation  of  the  form: 

dM^  du* 

ai  dS*  +  a2  dS*  =  a3  (19) 


E.D.  Kennedy,  Wake -Like  Solutions  to  the  Laminar  Boundary  Layer 
Equations,"  AIM  J.  2,  225-231  (1964). 
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where  a^,  and  a^  contain  no  S*  derivatives  (the  from  the  displace¬ 
ment  thickness  integrals  can  be  eliminated  by  using  the  momentum  equa- 
tion) .  Equation  (8)  has  the  same  form  (again  eliminating  using  ti 
momentum  equation) : 

dM  du* 

bld^  +  b2  d#-b3  <20) 

gM  du* 

e  c 

and  these  two  equations  solved  for  and  give: 


dM 

e 

a3b2 

‘a7b3 

h 

dS*  = 

alb2 

"  Vi 

D 

du* 

C 

*lb3 

"  a3bl 

_N2 

dS*  * 

alb2 

'  a2bl 

D 

Initially*  and  D  have  a  like  sign,  the  opposite  sign.  When 


is  too  large ^  changes  sign  first,  then  D  approaches  zero  so  that 

oo^  -oo.  When  M6o  is  too  small,  changes  sign  first  and 


the  centerline  velocity  decreases,  then  D  approaches  zero.  The  eigen 

solution  has  and  D  (and  therefore  N^)  changing  sign  simultaneously. 

It  should  be  noted  that  if  one  of e  equations  ^19)  or  (20)  were  to  be 

eliminated  by  specifying  either  (S*)  or  (S*) ,  the  singular 

behavior  would  no  longer  occur.  This  suggests  an  alternate  method  of 

solution  of  the  equations  diich  described  in  Appendix  I. 

Of  particular  interest  to  those  using  the  integral  methods  is  the 

profile  shape  obtained  by  the  finite  difference  method.  A  convenient 

way  of  presenting  the  shape  i-  to  plot  vS*  u>  as  Fig*  4#  The 

initial  profile  at  the  wake  stagnation  point  (u  0)  and  the  calculated 

profile  at  u  «  .3  u  are  shown.  The  initial  conditions  coincide,  with 
c  e 

one  of  the  Stewartson  family  (p  =  -.1988),  The  shapes  of  other  Stewart- 
son  profiles,  also  plotted  in  Fig.  4  (p  =  -.40,  -.30)  show  clearly  that 
this  family  very  accurately  represents  the  calculated  non-similar  pro- 
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files.  However,  it  should  be  remembered  that  the  use  of  a  Stevartson 
profile  shape  at  the  wake  stagnation  point  was  arbitrary  and  there  is 
some  doubt  that  such  a  choice  is  realistic. 

This  study  has  shown  that  this  particular  viscous  interaction 
problem  can  be  solved  in  a  straightforward  manner  using  a  finite  differ 
ence  method.  A  considerable  effort  has  been  made  to  extend  the  method 
of  calculation  to  regions  cr  reverse  flow,  without  success.  Several 
hybrid  methods  have  been  investigated,  in  which  integrals  of  the  conser 
vation  equations  are  used  in  the  reverse  flow  region  and  finite  d  Infer¬ 
ences  are  used  elsewhere.  All  of  these  attempts  failed  because  the 
profile  shapes  used  for  the  integral  method  did  not  have  the  proper 
qualitative  features.  This  sensitivity  to  the  details  of  the  profile 
shape  in  the  reverse  flow  region  was  also  pointed  out  by  Webb,  et  al^ , 
who  used  a  momentum  integral  method  over  the  entire  shear  layer.  The 
results  of  an  integral  method  are  in  doubt  whenever  such  sensitivity 
exists,  so  it  appears  that  additional  work  is  needed  in  the  region 
between  the  body  and  the  wake  stagnation  point  before  any  quantitative 
results  can  be  obtained  from  the  interaction  model. 
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APPENDIX  I 
INVERSE  METHOD 


Computing  through  a  saddle  point  singularity  is  difficult  even  if 
the  equations  are  simple  in  form  since  the  eigenvalue  can  be  Mund  only 
to  a  limited  number  of  significant  figures.  Assuming  that  one  is  then 
close  eno  igh  to  the  critical  point  to  make  some  kird  of  extrapolation 
valid,  the  calculation  can  be  restarted  on  the  other  side  and  proceed 
from  there  without  difficulty.  However,  when  a  partial  differential 
equation  is  included  in  the  set  of  equations,  the  computing  time 
involved  in  this  procedure  can  become  excessive.  A  considerable  saving 
in  time  would  result  fro^  a  method  which  did  not  require  that  such  a 
large  number  of  significant  figures  be  carried. 

The  method  which  will  be  descgi^ed  here  has  this  property,  and 
consists  essentially  of  guessing  anC*  so*v^nS  tbe  momentum 

equation,  using  Eq.  (8)  to  give  ,  ^The^  remaining  auxiliary 

equations  are  not  satisfied  (that  is  d( — —  )/dS*  from  Eq#  (11)  does 

P«Uoc£*  viSb  duc 

not  agree  with  A(  from  Eq,  (13),  and  a  new  guess  of  (S*) 

is  made  in  such  a  way  that  agreement  impioves  in  Succeeding  iterations. 

There  Ls  no  guarantee  of  convergence  and  the  key  to  the  method  is, 

obviously,  the  type  of  iteration  used. 

du* 

For  any  given  guess  of  (S*) ,  the  momentum  equat*  and  Eq,  (8) 

uo  dM 


can  be  solved  for  G(m, S*) ,  u*(S*)  and  ^ 


dS* 


(S*) .  The  initial  condi¬ 


tions  required  for  this  solution  are  G(Q,o)  and  u*(0)  .^The  values  of 

Me(S*)  do  not  appear  in  these  equations,  so  that  (S*)  is  found 

e 

without  specifying  anywhere.  The  initial  condition  needed  to  inte¬ 
grate  —  -jtt  (S*)  can  therefore  be  applied  at  any  value  of  S*.  The 

Mg  dhw 

eigen-solution  we  are  seeking  Ms  M  -  M  at  the  origin  and  M  -*  M 
&  a  e0  *  e  eco 

asymptotically  far  downstream.  Since  M  iCO  is  known,  we  will  integrate 

h - in  the  upstream  direction,  assuming  that  the  domain  has 

do* 

e 

been  chosen  large  enough  to  make  the  asymptotic  value  of  M£  reasonable 


at  the  downstream  end. 


Since  ^(co)ileaj  only  for  the  eigen-solution. 
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this  means  that  if  the  procedure  converges  at  all,  it  must  converge  to 
the  desired  solution. 

.  dM 

The  fact  that  the  differential  equation  and  —  -tXaT  (S*)  are 

Me  dS* 

integrated  in  opposite  directions  introduces  the  major  difficulty 

encountered  in  applying  the  method,  since  it  makes  it  impractical  to 

use  a  straightforward  iteration  scheme  such  as  Newton's  method.  If  they 

dug 

were  integrated  in  the  same  direction,  a  local  change  in  -  would  affect 

Eqs .  11  and  13  only  downstream  uf  the  point,  whereas,  when  they  are 

integrated  in  opposite  directions,  the  change  is  felt  downstream  due  to 

the  change  in  and  upstream  due  to  the  effect  on  all  preceding  values 

of  M  .  The  iteration  method  used  was  the  simplest  possible;  a  complete 
e  dug 

calculation  was  made  using  the  current  guess  of  — r(S*) ,  and  a  second  com- 

abs 

plete  calculation  was  made  using  1+s  times  the  current  guess  (where  €«1). 

The  results  of  these  two  calculations  were  used  to  make  a  point  oy  point 

dug 

linear  extrapolation  to  that  value  of  giving  agregmgnt  between 

equations  11  and  13.  All  of  the  indicated  changes  in  ttt  were  then 

aS* 

reduced  by  a  factor  which  decreased  the  largest  fractional  change  to 
P(P<1) a  There  is  no  mathematical  foundation  for  this  iteration  scheme, 
but  a  rigorous  extrapolation  is  so  much  more  difficult  that  it  was 
decided  to  try  this  first.  How  veil  the  method  actually  works  is  illus¬ 
trated  in  Figure  5.  The  problem  here  is  the  same  as  the  one  used  to 
demonstrate  the  direct  method  calculation,  where  the  eigenvalue  was 

found  to  be  =  6*55.  The  domain  was  taken  to  extend  from  S*=0  to 
e°  dug 

S*«1.5  (at  which  point  =6.0).  The  initial  guess  of  ^jx  was 

taken  to  be  0.3  everywhere.  Initially,  values  of  p  =  .2,  €  =  ,05  were 

used.  To  speed  convergence,  £  was  changed  to  .05  after  8  iterations. 

About  15  iterations  were  re 4  ired  to  find  the  eigen -value  to  the  same 

number  significant  figures  as  the  direct  method.  The  agreement  after 

Pc»U00  5* 

20  iterations  between  calculated  from  Eq,  11  with  that  calcu¬ 

lated  from  Eq^  13  is  shown  in  Fig.  6.  This  agreement  was  very  good 
near  the  origin,  and  was  poorest  at  about  S*=*9.  Further  iterations 
failed  to  improve  the  agreement  (which  is  hardly  surprising  considering 
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du*  dM 

the  iteration  method  used).  The  values  of  and  —  — f-  found  are 

dS*  M  dS* 

compared  with  the. direct  method  results  in  Figs.  7  anS  8.  The  agreement 

Pa»ucoO*  '  Pmum&* 

between  d(-j=-  )/dS*  from  quation  11  with  A(  )/AS*  from  Equation  13 

was  best  near  the  origin  and  was  poorest  at  about  S*=.7.  In  the  region 
between  the  origin  and  the  last  point  plotted  in  Figs.  7  and  8,  the  maxi¬ 
mum  discrepancy  between  the  two  values  is  10°/ o. 

It  is  interesting  to  note  that  the  direct  method  solution  could  not 
be  carried  beyond  S*=.3  when  four  significant  figures  were  required  of 
the  profiles,  while  the  indirect  method,  requiring  the  same  number  of 

significant  figures,  was  reasonably  accurate  to  S*  =  .5.  In  addition, 

PccUcofa* 

the  errors  in  agreement  between  d  (— )/dS*  and  the  finite  differences 
derivative  appear  to  be  distributed  in  such  a  way  that  the  displacement 
thicknesses  far  downstream,  when  calculated  in  the  two  different  ways  as 
shown  in  Fig.  6,  agree  quite  well. 

The  most  surprising  feature  of  this  method  is  that  it  works  at  all, 
even  without  considering  the  crudnessof  the  iteration  scheme.  The  biggest 
drawback  is  that  the  degree  of  convergence  is  limited,  but  it  is  possible 
that  this  can  be  improved  by  refining  the  iteration  procedure. 
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FIG.  1  BASF  FLOW  REGION 
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FIG.  4 

COMPARISON  BETWEEN  PROFILE  SHAPES 
FROM  THE  NON-S1HLLAR  CALCULATION 
AND  THOSE  FROM  THE  STEWARTSON- 
KENNEDY  SIMILAR  SOLUTION 


FIG.  5 
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FIG.  6 

CONVERGENCE  OF  INVERSE 
METHOD,  20  ITERATIONS 
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